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We describe and demonstrate an algorithm for the efficient simulation of generic stochastic dy¬ 
namics of classical degrees of freedom defined on the vertices of a locally tree-like graph. Networks 
with cycles are treated in the framework of the cavity method. Such models correspond for example 
to spin-glass systems, Boolean networks, neural networks, or other technological, biological, and 
social networks. Building upon ideas from quantum many-body theory, the algorithm is based on 
a matrix product approximation of the so-called edge messages - conditional probabilities of ver¬ 
tex variable trajectories. The matrix product edge messages (MPEM) are constructed recursively. 
Computation costs and precision can be tuned by controlling the matrix dimensions of the MPEM 
in truncations. In contrast to Monte Carlo simulations, the approach has a better error scaling 
and works for both, single instances as well as the thermodynamic limit. As we demonstrate at 
the example of Glauber dynamics, due to the absence of cancellation effects, observables with small 
expectation values can be evaluated reliably, allowing for the study of decay processes and temporal 
correlations. 

PACS numbers: 64.60.aq, 02.50.-r, 02.70.-c 


Introduction. - In the last years, we have seen increased 
efforts of statistical physicists to tackle stochastic dynam¬ 
ical processes in networks in order to study various phe¬ 
nomena [1, 2] such as ordering processes, the spreading 
of epidemics and opinions, synchronization, collective be¬ 
havior in social networks, stability under perturbations, 
or avalanche dynamics. 

A drastic simplification can be achieved when short cy¬ 
cles in the network, defined by interaction terms, are very 
rare. This is the case for locally tree-like graphs such as 
random regular graphs, Erdds-Reny graphs, and Gilbert 
graphs. For such random graphs with N vertices, al¬ 
most all cycles have length > log N such that their effect 
is negligible in the thermodynamic limit [3]. For static 
problems, this has been exploited in the so-called cavity 
method [4], where conditional nearest-neighbor proba¬ 
bilities are computed iteratively within the Bethe-Peierls 
approximation. The method was very successfully ap¬ 
plied to study for example equilibrium properties of spin 
glasses [4], computationally hard satisfiability problems 
[5, 6], and random matrix ensembles [7]. 

This big success has motivated the generalization of 
the cavity method to dynamical problems [8, 9], which is 
known as the dynamic cavity method or dynamic be¬ 
lief propagation. As the number of possible trajecto¬ 
ries and, hence, the computational complexity increase 
exponentially in time, applications have however been 
quite restricted to either very short times [8, 10], ori¬ 
ented graphs [8], or unidirectional dynamics with local 
absorbing states [9, 11-13]. In the latter case, one can 
exploit that vertex trajectories can be parametrized by 
a few switching times. Another idea has been to ne¬ 
glect temporal correlations completely as in the one-step 
method [8, 14-16] or to retain only some At = 1 corre¬ 


lations as in the 1-step Markov ansatz [17]. While this 
works well sometimes for stationary states at high tem¬ 
peratures, such approximations are usually quite severe 
for short to intermediate times or low temperatures. 

In this paper, we present an efficient novel algorithm 
for the solution of the parallel dynamic cavity equations 
for generic (locally tree-like) graphs and generic bidirec¬ 
tional dynamics. The central objects in the dynamic cav¬ 
ity method are conditional probabilities for vertex trajec¬ 
tories of nearest neighbors - the so-called edge messages. 
As temporal correlations are decaying in time and/or 
time difference \t — t'\, we can approximate each edge 
message by a matrix product, i.e., there is one matrix 
for every edge, edge state, and time step, encoding the 
temporal correlations in the corresponding part of the 
evolution. It turns out that the dimensions of these ma¬ 
trices do not have to be increased exponentially in time. 
One can obtain quasi-exact results with much smaller 
matrix dimensions. Computation costs and precision can 
be tuned by controlling the dimensions in truncations. 
The idea of exploiting the decay of temporal correla¬ 
tions to approximate edge messages in matrix product 
form is in analogy with the use of matrix product states 
[18-22] for the simulation of strongly correlated, mostly 
one-dimensional, quantum many-body systems. These 
have been used very successfully in algorithms like the 
density-matrix renormalization group [23, 24] to compute 
for example quantum ground-state properties, often with 
machine precision [25]. Besides lifting the restrictions of 
the aforementioned approaches, the matrix product edge- 
message (MPEM) algorithm can also outperform Monte 
Carlo simulations (MC) of the dynamics in important re¬ 
spects. In particular, besides allowing for the simulation 
of single instances, alternatively, one can work directly in 
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the thermodynamic limit. Perhaps more importantly, it 
has a favorable error scaling. While statistical errors in 
MC decay very slowly with the number of samples Ns as 
MPEM yields also observables with absolutely 
small expectation values with very good precision which 
is essential for the study of decay processes and temporal 
correlations. 

The dynamic cavity method. - Let denote the state 
of vertex i at time step t, and cr* := (cr*, cr^i ■ • ■) the 
state of the full system at time t. Given the state prob¬ 
abilities P{cr*) for time t, we evolve to the next time 
step, P(cr*+i) = W{cr*'^^\cr*)P{(T*), by applying the 

stochastic transition matrix W. As vertex i interacts only 
with its nearest neighbors j € di, the probability for 
only depends on the states crj of these vertices at the pre¬ 
vious time step such that the global transition matrix W 
is a product of local transition matrices Wi, 

(1) 

i 

Here = 1) cr^- is the state of the 

nearest neighbors of vertex i at time t. In the cavity 
method [4, 8, 9], one neglects cycles of the (locally tree¬ 
like) graph according to the Bethe-Peierls approximation 
to reduce this computationally complex evolution to the 
dynamic cavity equation [8, 9] 
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which only involves the so-called edge messages /r for the 
edges of a single vertex i. For simplicity, we have assumed 
that vertices are uncorrelated in the initial state such that 
The edge messages in 

the dynamic cavity equation (2) are conditional probabil¬ 
ities for the trajectories tr* := (tT°, (7,^,..., cr|) and 
on edge (i,j). Specifically, if we consider a tree graph 
and cut off everything “right” of vertex j as indicated 
in Figure 1 by the dashed line, denotes 

the conditional probability of a trajectory ct* on vertex i, 
given the trajectory on vertex j. Eq. (2) constructs 
out of the edge messages 

of the previous time step. This is exact for tree graphs 
and covers locally tree-like graphs in the Bethe-Peierls 



Figure 1: (Color online) Part of a locally tree-like interaction 
graph with vertex degrees z = 3. 


approximation. Although we have gained a lot in the 
sense that the computational complexity is now linear in 
the system size, it is still exponential in time t, if we were 
to encode the edge messages without any approximation. 

Canonical from of an MPEM. - To circumvent this 
exponential increase of computation costs, we can exploit 
the decay of temporal correlations and approximate the 
exact edge message by a matrix product 


t-i 

1) = A®^.(a°) [ n dl(i(ar V|) 

X (3) 

The particular choice of assigning vertex variables {cr®} 
and {cr®} to the Ms x matrices A-^^ ((T®“^|cr®) oc¬ 

curring in the matrix product (3), is advantageous for the 
implementation of the recursion relation (2) for MPEMs 
as will become clear in the following. In order for the 
matrix product to yield a scalar, we set Mq = Mt +2 = 1- 
MPEM evolution. - The time-evolution starts at t = 0 
with = pi{a^). Using the dynamic cavity equa¬ 

tion (2), we iteratively build matrix product approxima¬ 
tions for edge messages for time t -I- 1 from those for time 
t. It is simple to insert the matrix product ansatz (3) 
for the edge messages in the dynamic cavity equation, 
but not trivial to bring the resulting edge message again 
into the canonical MPEM form as required for the sub¬ 
sequent evolution step. The specific assignment of the 
vertex variables to matrices in Eq. (3) has been chosen 
such that all contractions (products and sums over ver¬ 
tex variables) occurring in the cavity equation are time- 
local in the sense that, given MPEMs in 

canonical form for all neighbors k € di\{J}, the resulting 
/ii_>j(d’^~''^|d*) can be written in (non-canonical) matrix 
product form as 




l-p= 


(0) 

i^j 


. i+1 




) . ( 4 ) 


for 1 < s < t 

are obtained by contracting a local transition matrix 


As depicted in Figure 2b, the tensors 




_s-l> 


with tensors from the time-t MPEMs. 


This contraction entails a sum over the z — 1 common 
indices where z = \di\ is the vertex degree. As¬ 

suming for the simplicity of notation that the matrix di¬ 
mensions Ms for all time-t MPEMs are identical, the re¬ 
sulting matrices (cr®|cr®“^) 


0 


kGdi\{j} 


di\{j} 

have left and right indices 


\z-l 


of dimensions Ms = {MsY ^ and Ms+i = {Ms+i) 
respectively. The contraction for tensor (7^+^) is very 
similar and = pi{a°) ■ 

MPEM truncation. - In preparation for the next time 
step, we now need to bring the evolved edge message (4) 
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Figure 2: (Color online) (a) Graphical representation of a matrix product edge message in canonical form (3). Connecting 
lines indicate summations over indices, (b) For each time step (2), tensors of the evolved matrix product in 

Eq. (4) are built by contracting the local transition matrix Wi with MPEM tensors of messages fik^i, incident to vertex i, 
where k £ di \ {j} = {fci,..., and a := (ai,..., Oz-i). (c) Evaluation of probabilities P{crj,aj) as in Eq. (8). 


back into canonical form (3). Furthermore, we need to 
introduce a controlled approximation that reduces the 
matrix dimensions because they would otherwise grow 
exponentially in time. Both, a reordering of the vertex 
variables (t| and cr® in the matrix product (4) and a con¬ 
trolled truncation of matrix dimensions can be achieved 
by sweeping through the matrix product and doing cer¬ 
tain singular value decompositions (SVD) [26] of the ten¬ 
sors. 

Let us shortly explain the notion of truncations at the 
example of a matrix product 7 (n) := Aq° , 

where A^= is an Mg x matrix and Mq = Mt+i = 1. 

In order to reduce in a controlled way, e.g., the left matrix 
dimension of we suggest to do an SVD of the 
matrix product such that 

7(n) =: TriL,nii — '^^Yn^^k^kZk,rLii (5) 

fe=l 

where we have grouped the variables to := 

(no,..., riu-i) and := (n„,..., rit). Y and Z are iso¬ 
metric matrices; Y^Y = 1 and = 1. Now, truncat¬ 
ing some of the singular values Ai > A 2 > • • • > Am„ > 0, 
such that only the M'^ largest are retained, we obtain the 
controlled approximation 

7trunc(^) •— ^ ^ Y^j^ k^kZk^nii 

k<M^ 

with error ||7 - 7trunc||^ = ^ A^. (6) 

k>M'^ 

Note that this truncation scheme guarantees the mini- 

1 /2 

mum possible two-norm loss jlAyH = 
for the given new matrix dimension AI'^. 

While it seems very desirable to discard unimportant 
information and control the growth of computation cost 
through such truncations, the SVD (5) appears to be an 
insurmountable task. Assuming that each variable rig can 
take d different values and that 2u < t+1, the cost for the 
SVD would scale exponentially in time like [26]. 

However, the beauty of matrix products is that such an 


SVD can in fact be done sequentially with linear costs of 
order tdM^ {M := maxg Mg) as follows. First, we do an 
exact transformation of the matrix product to bring it to 
the form 7 (n) = Y^° ■ ■ ■ • • • Z”*, where 

tensors Yg and Zg obey the left and right orthonormality 
constraints 

^(Vg7n" = l and ^Yg"(Yg7 = l, (7) 

n n 

respectively. This is achieved through a sequence of 
SVDs. It starts with the SVD =: Y^^AoVo, 

where Ag is a diagonal matrix of singular values, Vg 
is isometric according to VqV^ = 1 and Vg obeys 
Eq. (7). The sweep continues with the SVD AgVgA"^ =: 
Y^^AiVi and so on until the computation of Vx-i- 
Analogously, we do a second sequence of SVDs start¬ 
ing from the right and ending with A^'^^Uu+ 2 Au +2 ='■ 
Uu+iAu+iZJ^^ and, finally, define the central tensor 
as := A„_it4_iA]];"[/u+iA„+i. After this some¬ 
what laborious preparation, we can do the actual trunca¬ 
tion, based on the SVD =: with the same 

singular values Ai > ••• > Am„ as in Eq. (5). With 
the Mu X M'u matrix [Atruncjfefc' := 4fc'Afc, the trun¬ 
cated matrix product (6) takes the form 7trunc('>T') = 

-^"11-2 /TZ^u-l 7-T A \ yriu 

Jg '''^u-2 V^tt-1 ^uMruncj'^u ' " ■ 

With this tool in hand, we can now truncate the 
evolved edge message (4) and bring it back into canon¬ 
ical form (3). In a first sweep from right (s = t -I- I) 
to left (s = 0), using SVDs, we can sequentially im¬ 
pose the right orthonormality constraints [see Eq. (7)] 
on the C-tensors. In a subsequent sweep from left 
to right, again based on SVDs, we can now truncate 
the tensors to decrease bond dimensions from Mg to 
something smaller. What is left, is to reorder the in¬ 
dices {ct|} and {(t|} of the vertex variables. In a 
sweep from right to left, we go from the variable assign¬ 
ment io'i){crl\aj)... (cr-'''^|cr*) in the truncated and or- 

thonormalized version Illii of 

the MPEM (4) to the assignment (77) ■ ■ • (7l7))^i^^) 
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Figure 3: (Color online) (a) Magnetization and, (b), connected temporal correlations for Glauber dynamics on 2 = 3 random 
regular graphs of different sizes for MC and in the thermodynamic limit for the MPEM and 1-step Markov approaches. Because 
of odd-even effects in the dynamics, only even time steps are shown. For MC (with Ns samples), the errors of the magnetization 
[lower panels in (a)] are quantified by the standard deviation of the magnetization, i.e., under ignorance of remaining finite-size 
effects. For MPEM and 1-step Markov, errors are quantified by the deviation from the result of the most precise (quasi-exact) 
MPEM simulation (truncation threshold Atrunc = 10~® for /3 = 1 and Atmnc = 10“^ for )? = 1/4). In plot (b), the three MPEM 
curves for Atrunc = lO”'^, 10“®, 10“® overlap up to time t = 24. 


in the matrix product 

s^O 

At the right boundary, we start with an SVD and 
controlled truncation (ct‘+V j) 

and continue with C|!|j-(cr*|cry^) 

xC/(‘+i)(aj)y‘y) «: C/(‘)(ayi)yL,i4l!l^.(a‘|ap and 
so on until ending at s = 0. In an analogous final 
sweep from left to right, we change to the canonical vari¬ 
able assignment (cr°)(cr°|crj)... as in 

Eq. (3). After executing these steps for all edge messages, 
the next evolution step from < -|- 1 to t -I- 2 can follow. 

Evaluation of observables. - The joint probability of 
trajectories aj and aj~^ for the vertices of an edge {i,j) 
is given by the product of the two corresponding edge 
messages. After marginalization, one obtains for example 
the probability for the edge state ((T|,cr*) at time t as 

P(a‘,a*)= ^ Mw,(d‘|dyi)/r,^.(d‘|d*-i). (8) 

-t-1 -t-1 

In the MPEM approach, this can be evaluated efficiently, 
as indicated in Figure 2c, by executing the contractions 
sequentially from left (s = 0) to right {s = t — 1). Simi¬ 
larly, one can for example also compute temporal corre¬ 
lators {crlaf) from probabilities P{crl,af). 

Exemplary applieation. - Figure 3 compares the sim¬ 
ulation of Glauber dynamics [27] using our MPEM al¬ 
gorithm to MC simulations and to the 1-step Markov 


approximation [17]. Specifically, we have Ising spins in¬ 
teracting ferromagnetically on 2 = 3 random regular 
graphs, with local transition matrices = 

exp(/3 CTy^CTp/.Z. In the initial state, all spins 

have magnetization (cr°) = 1/2, i.e., Pi{f) = 3/4. Be¬ 
sides being applicable for single instances of finite graphs, 
the MPEM approach gives also direct access to the ther¬ 
modynamic limit. For disordered systems this can be 
done in a population dynamics scheme. The homoge¬ 
neous case, considered here, is particularly simple as all 
edges of the graph are equivalent in the thermodynamic 
limit. Hence, one can work with a single MPEM. Fig¬ 
ure 3a shows the evolution of the magnetization. In 
the ferromagnetic phase (/3 = 1), it approaches a fi¬ 
nite equilibrium value, whereas it decays to zero in the 
paramagnetic phase {j3 = 1/4). As shown for /? = 1, 
MC simulations contain finite size effects which become 
small for the system with 2048 sites. MC errors decrease 
slowly when increasing the number of samples as 
This is problematic for observables with small 
absolute values where cancellation effects make it diffi¬ 
cult to get a precise estimate. This is, e.g., apparent 
in the magnetization decay for /3 = 1/4 which, in con¬ 
trast, is very precisely captured with MPEM. In these 
simulations, we control the MPEM precision by keeping 
only singular values A^ above a threshold, specified by 
^k/J2k' > Atrunc- Decreasing Atmnc, increases preci¬ 

sion and computation costs. The 1-step Markov approxi¬ 
mation [17] is not suited to handle temporal correlations. 
At long times it performs well for /3 = 1/4 and fairly 
good for /3 = 1, but deviates rather strongly at earlier 
times. Figure 3a shows the connected temporal correla- 
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tion function — {crl){crl) for /3 = 1 as a function of 

t — s for several times. Its decay behavior can be difficult 
to impossible to capture with MC. In the example, MC 
deviations are often orders of magnitude above those of 
the numerically cheaper MPEM simulations. 

Conclusion. - The novel MPEM algorithm, based on 
matrix product approximations of edge messages allows 
for an efficient and precise solution of the dynamic cav¬ 
ity equations. Besides lifting restrictions of earlier ap¬ 
proaches, mentioned in the introduction, it gives direct 
access to the thermodynamic limit, and its error scaling 
is favorable to that of MC simulations. We think that it 
is a very valuable tool, particularly, as it yields tempo¬ 
ral correlations and other decaying observables with un¬ 
precedented precision and gives access to low-probability 
events. This opens a new door for the study of diverse 
dynamic processes and inference or dynamic optimiza¬ 
tion problems for physical, technological, biological, and 
social networks. 

We thank G. Del Ferraro for providing data to bench¬ 
mark our 1-step Markov simulation for Figure 3 and ac¬ 
knowledge supported by the Marie Curie Training Net¬ 
work NETADIS (FP7, grant 290038). 
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